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Abstract 

We study quench dynamics of the Bose-Hubbard model by exact diagonahzation. Initially the 
system is at thermal equilibrium and of a finite temperature. The system is then quenched by 
changing the on-site interaction strength U suddenly. Both the single-quench and double-quench 
scenarios are considered. In the former case, the time-averaged density matrix and the real-time 
evolution are investigated. It is found that though the system thermalizes only in a very narrow 
range of the quenched value of U, it does equilibrate or relax well in a much larger range. Most 
importantly, it is proven that this is guaranteed for some typical observables in the thermodynamic 
limit. In order to test whether it is possible to distinguish the unitarily evolving density matrix from 
the time-averaged (thus time-independent), fully decoherenced density matrix, a second quench is 
considered. It turns out that the answer is affirmative or negative according to the intermediate 
value of U is zero or not. 

PACS numbers: 05.70.Ln, 05.30.Jp, 05.30.-d 
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I. INTRODUCTION 



Out-of-equilibrium dynamics following a quantum quench is a topic of intense study 
at present. The theme is pursued primarily along two lines. The first one is about the 
equilibration and thermalization mechanism of a quantum system iMlOj]. a fundamental yet 



still open issue in statistical physics. The second one is about the the real-time dynamical 



behavior of a many-body system 
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iMlSl]. which is highly non-trivial in the regime where the 
quasi-particle picture breaks down. 

Among all the models investigated so far, the Bose-Hubbard model takes a special posi- 
tion. As a paradigmatic strongly-correlated model, it can be realized accurately with cold 
atoms in optical lattices, and especially, the parameters can be controlled (e.g. changed 
suddenly) to a high degree 16l-ll8|. This nice property makes it an ideal candidate to in- 
vestigate quantum quench dynamics both theoretically and experimentally. Up to now, in 



13|, 



the few theoretical works on the quench dynamics of the Bose-Hubbard model J 5, 
the state of the system before the quench is always assumed to be the ground state of the 
initial Hamiltonian. That is, the system is assumed to be at zero temperature initially. How- 
ever, in this paper we shall start from a thermal equilibrium state. One should note that 
this scenario is actually more experimentally relevant. Because in current experiments, one 

generally gets not a single tube of cold atoms, but instead a two-dimensional array of one- 

j 

dimensional lattices for the cold atoms |18|. In other words, an ensemble of one-dimensional 
Bose-Hubbard models is obtained in one shot. Moreover, in view of the fact that the cold 
atoms are at finite temperatures necessarily 
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20|, it is reasonable to start from a thermal 
state described by a canonical ensemble density matrix [see Eq. below]. 



As emphasized by Linden et al. 



2lj, in the pursuit of thermalization, it is important 



to distinguish the two closely related but inequivalent concepts of equilibration and ther- 
malization. The latter is much stronger and has the trademark feature of the Boltzmann 
distribution, whereas the former refers only to the stationary property of the density matrix 
of a (sub)system or some physical observables. It is highly possible that a system equilibrates 

or the Bose-Hubbard model. As re- 



5| and in the present paper (finite 



but without thermalization. This is actually the case 
vealed both in previous works (zero temperature case) 
temperature case), the Bose-Hubbard model thermalizes only if the quench amplitude is not 
so large, at least at the finite sizes currently accessible. However, it will be shown below 
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that in a much wider range of parameters, some generic physical observables equihbrate very 
well. Among them are the populations on the Bloch states, which are ready to measure by 



the typical time-of- flight experiment j22| . Remarkably, this is actually guaranteed for these 
quantities in the thermodynamic limit, i.e., when the size of the system gets large enough. 

The equilibration behavior of the physical observables imposes a question. It is ready to 
recognize that the equilibration of the physical observables is largely an effect of interference 
cancelation. It never means that the density matrix has suffered any dephasing or decoher- 
ence. Actually, the density matrix evolves unitarily and in the diagonal representation of the 
Hamiltonian, its elements simply rotate at constant angular velocities. A natural question 
is then, does the time- dependence of the density matrix has any chance to exhibit it, given 
that it is almost absent in the average values of the physical observables? This leads us to 
consider giving the system a second quench. The concern is, would the system yield different 
long-time behaviors if the second quench comes at different times? It turns out that the 
answer depends on whether the intermediate Hamiltonian is integrable or non-integrable. 
In the former case, the density matrix shows repeated appreciable recurrences and thus the 
dependence on the second quench time is apparent. In the latter case, on the contrary, the 
density matrix shows no sign of recurrence and quantitatively similar long-time dynamics is 
observed for quenches at different times. 

This paper is organized as follows. In Sec. [Tll the setting of the problem and the basic 
approaches are given. In Sec. IIIIl the dynamics after a single quench is studied. The 
time-averaged density matrix and the real-time evolution of some physical observables are 
investigated in detail. Based on the observation in this Section, we proceed to study the 
scenario of a second quench in Sec. IIVI Finally, we summarize the results in Sec. |Vl 

II. BASIC FORMALISM 

The time-dependent Hamiltonian of the Bose-Hubbard model is {h = ks = 1 throughout 
this paper) 

H{t) = - J ^(afa^+i + a\^^ai) + — ^ ^ a]alaiai. (1) 
1=1 1=1 

Here M is the number of sites (the total atom number will be denoted as A^) and aj (a;) is 

the creation (annihilation) operator for an atom at site /. Note that here periodic boundary 
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condition is assumed. The parameters J and U are the nearest-neighbor hopping strength 
and the on-site atom-atom interaction strength, respectively. Note that the dynamics of 
the system depends only on the ratio UjJ, thus we will set J = 1 throughout. We say 
the system is quenched if U is changed suddenly at some time from one value to another 
value. Experimentally, for cold atoms in an optical lattice, this can be realized by using the 
Feshbach resonance. 

Assume that initially the parameter U is of value Ui (the corresponding Hamiltonian 
is denoted as Hi), and the system is at thermal equilibrium and of inverse temperature 
Pi = 1/Tj. Denote the m-th eigenvalue and eigenstate of Hi as E^^ and \ipln}j respectively. 
The initial density matrix of the system is then 

1 ^ 
p. = -exp(-Ai^.) = (2) 

' m=l 



where = ^^=1 ^^P(~/^«-^m,) partition function and = exp(— is the 

probability of occupying the eigenstate IV^J^). Note that D = is the dimension of 

the Hilbert space H. The density matrix at time t is given formally as p{t) = U{t)piW{t), 
with U{t) = Texp[-i jldrH^T)]. Here T means time ordering. 

The Hamiltonian H{t) is invariant under the translation (a/,a|) — )■ (a/+i,a|^^). This 
indicates that the total quasi-momentum of the system q = Xlfclo^ ka\ak (mod M), where 
a\ = ^i^i exp{i2Trkl / M)al is the creation operator for an atom in the A;-th Bloch state, 
is conserved. This property implies that if the full Hilbert space is decomposed into M 
subspaces according to the values of q, i.e., "H = (BgiiQ^'H^''\ the Hamiltonian and the density 
matrix are always block-diagonal with respect to the g-subspaces, i.e.,H{t) = (Bg^^H^'^^t) 
and p{t) = ©^lo^p'''^''(t) 23|, l24l- It is then possible to study the dynamics in each subspace 
individually (which saves a lot of computational resource) and then gather the information 
together (note that for the expectation values of quantities like a|,afc, there are contributions 
from each subspace). Here it is necessary to mention that though we should have done 
the gathering or averaging process for many quantities studied below, we would rather not 
do so, because it is observed that the system behaves quantitatively similar in all the q- 
subspaces {25 1. A single g-subspace captures the overall behavior very well. Therefore, our 
strategy is to focus on some specific g-subspace (g = 1 actually) and take the normalization 
tr{p^'^\t)) = 1. It is understood that in the following all Hamiltonians, density matrices, 
eigenvalues, and eigenstates refer to those belonging to this specific q-subspace. We will drop 
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the superscript q for notational simplicity. 



III. A SINGLE QUENCH 

Suppose at time t = the system is quenched by changing the value of U from Ui to 
f//^, which is then held on forever. The Hamiltonian later will be denoted as i^/j, and the 
eigenvalues and eigenstates associated will be denoted as E^^ and \ipn), respectivley. In the 
representation of {lipl^)}, the density matrix at time t is then simply (in this paper (■ ■ ■) 
means quantum state averaging while means time averaging) 

Pit) = 5^ e-^(^--^"^)*(^^^|p.|^„^^)|^4^)(^;^^|, (3) 

m,n=l 

where Dg D/M is the dimension of the specific g-subspace. It will prove useful to define 
the time-averaged density matrix 

1 '■^ 



Da 



E {i^!;i\pM'M!;i){^i'\- (4) 



m,n=l 
Tpli — pi'i 



The time- averaged density matrix is of great relevance for our purposes. First, it is both 
time-independent and variable-independent. Second, the time-averaged value of an arbitrary 
operator O is given simply by (O) = limr__i.oo Jq tr{p{t)0)dt = tr{pO). That is, the time- 
averaged density matrix contains the overall information of the dynamics of the system. 
Actually, as we will see later, for some quantities which fluctuate little in time, the time- 
averaged density matrix tells almost a complete story. Third, the process of averaging over 
time is a process of relaxation in the sense that the entropy associated with p is definitely 
no less than that with the density matrix at an arbitrary time, i.e., S{p) > S{p{t)) = S{pi). 



This is a corollary of the Klein inequality 



and is reasonable since pi contains all the 



information of p while the inverse is invalid. The equality also means that p{t) will never be 
damped, and time-averaging is essential. 

Note that when Uf^ ^ 0, generally there is no degeneracy between the eigenvalues of if/j. 
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Therefore the time-averaged density matrix is simply diagonal in the basis of {{ipn)}, i.e., 

m=l 

with 

1 

n=l 

being the population on the eigenstate l'?/'^). In the special case of Uf^ = 0, the Hamiltonian 
reduces to Hf-^ = J2k^=o^ with cok = —2Jcos{27ik/M). In this case, each eigenvalue 

is of the form '^/^rikUJk, under the constraints Xlfc'^fc ~ ^ Ylk = Q (mod M), and 
there can be level degeneracy. However, we can always make some unitary transforms in 
each degenerate subspace to make sure that p is in the form of ([5]). 



A. Time-averaged density matrix 

Since the time-averaged density matrix provides an overall information of the dynamics 
of the system, we look into it first. In Fig. [H we consider the scenario of starting from the 



same initial condition {Ui = 1, /3j = 0.3) but quenching to six different values of Uf^ (26 1. 
In each panel, the logarithms of Pm are plotted against the eigenvalues (red dots). We 
have compared p with a canonical ensemble density matrix pc, which is defined as 

^ tr(e-^/i^/i) 

under the condition tr{pcHf^) = tr{pHf^) = tr{piHf-^). Here the final inverse temper- 
ature, is the only fitting parameter. In Fig. [H the green dots which form a straight line 
correspond to Pc- 

We see that p exhibits many interesting features. In the case of f/j^ = 0, p agrees well 
with Pc throughout the spectrum. In the case of f//^ = 2, p agrees well with pc in the lower 
part of the spectrum, while deviates from it significantly in the higher part of the spectrum. 
But overall the two are in good agreement since the weight of the higher part is small. The 
case oiUf^ = —1 is somewhat the reverse of the Uf^ = 2 case. It is in the lower part of the 
spectrum that Inp^ fluctuates wildly. In the higher part Inp^ goes almost linearly. Since 
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(a) U,, =-10, =-(l.(]142,-) 




(d) U,, =2. =0.24343 





(c) r/j, =0, ,% =0.30802 




(c) C/,, =5, ,% .0.11653 



-15 -10 -5 



(f) Uf, =10, =0.048254 




FIG. 1: (Color online) Semilog plots of Pm versus the eigenvalues Em (red dots). The initial state 
is the same for all the figures, with parameters (M, A^, q, Dg) = (9, 9, 1, 2700), Ui = 1, and /3j = 0.3. 
The quenched values of U and the fitting inverse temperatures /3/^ are shown in the inserts. For 
comparison, the data with pc (green dots) and p'^ (blue dots) are also shown. The black lines at 
the bottom depict the coarse-grained density of states of Hf-^ (just for reference, not corresponding 
to the vertical axis). 



the weight is dominated by the lower part, pc is not a good approximation of p. In the 
strong interaction limits of t//^ = ±10, another feature takes the place. As a whole the red 
dots do not fall close to a single straight line, but they do form some stripes, and the stripes 
are almost parallel with a common slope close to (3i. It is easy to recognize that each stripe 
corresponds to a bump in the density of states of Hf-^ . 

In order to understand the various features in Fig. [1], we rewrite pm as 



Pm y 



(8) 



where Pm{E) = J2n l(V'nlV'm)P'^(-^ ~ ^n) ^ probability distribution [27| associated with 
\ipm)- Note that Pm{E) is an intrinsic property of independent of We have tried 
to characterize the distribution Pm{E) by its mean /x^, its second central moment cr^, and 
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FIG. 2: (Color online) The parameters /Xm, and Km [see Eq. ([9])] characterizing the probability 
distributions PmiE) associated with the eigenstates of Hf^. Note that Figs. [2^-[2]f correspond to 
Figs. [T^-[l]f, respectively. 



its third central moment k^, which are defined as follows, 



dEPUE) = {iJ^^m^pt), 

J dEPm{E){E - l^r, 
J dEP^{E){E - fir, 



^■m ) ) 



(9a) 
(9b) 
(9c) 



These quantities are presented in Fig. [2J These data enable us to understand Fig. [H Suppose 
for a distribution Pm{E) with {^m,o'm), we define a Gaussian distribution 



PLiE) 



exp 



{E-fir, 



(10) 



which shares the same mean and variance with P^ but has vanishing third central moment. 
Replacing Pm in Eq. (|8]) by P^, we get an approximation of pm, 

1 / „ 1 



P'm = Y. {-Pif^rn + 2 A^^m 



(11) 



In Fig. [H p'r^ are represented by the blue dots. We see that as a whole p'^ is a good approx- 
imation of Pm, except at the lower part of the spectrum in Fig. [T)d. The reason is clear — the 
Km's there are the largest throughout all the figures, which indicates that the corresponding 
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distributions Pm are wide and asymmetric and thus cannot be well approximated with a 
Gaussian distribution. 

Now we can understand the good fittings in Figs. ^ and [T]i. In these two cases, l^rn 
almost a linear function of E^, and cr^ does not vary so much, therefore the exponent in 
Eq. (ITT]) goes almost linearly with E^. The situation is similar in the higher part of the 
spectrum in Fig. [2)d, and therefore we have a good linear fitting for the higher spectrum 
part in Fig. In contrast, in Fig. [2^, /i^ varies wildly for adjacent E^, therefore we see in 
Fig. [1^ large fluctuations about the straight line. As for the parallel stripes in Figs. [1^ and 
[T]F, they are also understandable in terms of Figs. [2^ and [2f, where fim form parallel stripes. 
It is numerically checked and can be argued that the slopes of the stripes are almost unity. 
Actually we have 

E^^ = i^tlHfMt) 

= i^ljl^m^ljf^) + {Uf, - U.){^f^\H,ntW^), (12) 

where Hint = | Xlfii '^I'^!'^''^'- Note that in the limit of large \Uf^/J\, the kinetic term in 
the Hamiltonian ([1]) can be viewed as a perturbation to the second interaction term. The 
spectrum of the latter is highly degenerate and consists of integral multipliers of Uf^. The 
effect of the perturbation is to mix up the eigenstates of the interaction Hamiltonian with 
different eigenvalues and smooth the spectrum. That is why there are bumps in the density 
of states in Figs. [1^ and [1]F and two adjacent bumps are placed roughly Uj^ apart. By 
perturbation theory, it is easy to show that the second term in Eq. (fT2|) varies on the order 
of J^/|f//J ^ J among eigenstates belonging to the same bump. Therefore, approximately 
we have /x^ = Em " const for each bump and this explains why the stripes in Figs. [2^ and 
[2]F are of slope unity. In turn it explains [with the help of Eq. ( ITTl) ] why we have the parallel 
stripes in Figs. [1^ and[l]F, and especially the slopes are approximately 

It seems in Fig. [1] that pc is a good approximation of p only when \Uf^ — Ui\ is small. 
In Fig. [31 we employ the tools of distance D, fidelity F, and relative entropy Srei (for the 



definitions see [28|) between two density matrices to quantify the difference or resemblance 
between pc and p. There it is clear that only in the range of \Uf^ — Ui\ < 1, we have 
(D, 1 — -F, Srei) ^ 1, which means p is close to pc- In the subsequent subsection we will see 
that only in this range the expectation values of some generic physical observables according 
to p and Pc agree well. 



FIG. 3: The distance D and fidelity F between pc and p, and the relative entropy of pc with respect 
to p, as functions of Uj-^ . The initial state is the same as in Fig. [TJ 




5 10 15 20 25 30 35 40 20 40 60 80 100 120 5 10 15 20 25 30 35 40 

t t t 



FIG. 4: (Color online) Time evolution of the populations on the Bloch states {a\ak). The figures 
correspond to those in Fig. [1] in a one-to-one manner. In each figure, from up to down, the five 
lines correspond to = 0, . . . , 4. Another fc's are now shown because {a\ak) and {a\.^_^aM-k) are 
close to each other all the time. For each line, the markers of the same color on the right hand 
side indicate the average value predicted by p (*) or value predicted by pc (□), respectively. Note 
that in (b) and (e), the time span investigated is longer than that in others. This is because the 
transient times in (b) and (e) are relatively longer. 

B. Time evolution 

We now proceed to study the time evolution of the system after the quench. In Fig. [2l 
we show the time evolution of the populations on the Bloch states {a\ak)- The six sub- 
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figures correspond to those in Fig. [T] respectively. For all the f//i's and all the /c's, {(i\ak) 
equilibrate to their average values after a transient time, which is relatively longer in the 
cases of f//i = — 1 and 5. In the special case of Uf^ = 0, there is no fluctuation at all. The 
reason is simply that in this case, a^ak are conserved. We see that the time-averaged values 
of (a\ak) predicted by p (*) and pc (□) agree relatively well in the cases of f//^ =0 and 2. 
This is consistent with the closeness between p and pc for these two values of Uf^ , as revealed 
in Fig. [1] and Fig. [31 Here we would say the system thermalizes well in the Uf^=2 case, 
however, we would refrain making the same statement for the t//^ = case. The reason will 
be clear in the next Section. 

Figure mis about a finite-sized system with some specific initial condition. However, here 
we have some general statements. We argue that in the thermodynamic limit (M, N ^ oo 
with N/M fixed), as long as initially the system is at finite-temperature thermal equilibrium 
and described by a canonical ensemble density matrix as ([2]) , we should see steady behaviors 
of the physical variables like a^Ofe. 

Let A = a\ak and let A = J2mn^rnn\ipm){i^t\ ^^e representation of The 
ensemble-averaged value of A at time t is 

("(t) = P^^nAnm exp[-t{E^ - El')t], (13) 

mn 

where pmn = {i^mlPili^t) ■ Its time-averaged value is 

a = Urn — dta{t) = y^PmmAmm- (14) 

m 

Here note that for a generic Hamiltonian Hf^ , there is no level degeneracy. The time-averaged 
value of a^{t) is ^291] 

a? = lim — / dta^it) 

^ ^ PmmAmmPppApp -\- ^ ^ pmnAnmPnmAmn 
mp m^n 

^ ^ PmmAmm ^ ^ PppApp -\- ^ ^ |Pmn| \An 



'-mn 



p m^n 

= (I -\- ^ ^ I Pmn I l^mnl • (l^) 
m^n 

Note that here it is assumed that there is no degeneracy of energy gaps. Thus we have for 
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the variance of a{t) in time, A^a = — a^, 

A O ^ ^ |Pmn| l^mnl ^ ^ ^ |Pr?ira| |^mn| • (-^^) 
m^n mn 

Since A is semi-positive definite and bounded, we have |v4mnP < ^mm^nn < N'^- Thus we 
have 

mn 

Here we note that the summation is the square of the Frobenius norm of pi in the represen- 
tation of {iV'm)}) which is invariant in all representations and is preserved by an arbitrary 
unitary evolution |30|. Explicitly, we have 

mn mn 

m m 

= Y^^I^mf- (18) 

m 

We argue that this quantity, which depends only on the initial state, decays exponentially 
with the size M. Let E^^ increase with m. We have 

Q-PiE\ ^-PiE{ g-ftoiM 
m ' 

as M — )■ oo. Here in the ~ relation we used the fact the ground state energy E\ of Hi scales 
linearly with M and so does the free energy Fi of the initial state 3l|]- The coefficients a 
and 7 are independent of M. Moreover, it is easy to see that a > 7 for any with the 
equality taken only in the limit of /3j = +00 or Tj = 0^, and a; — 7 increases monotonically 
with Tj. This makes sure that p\ would not grow exponentially with M and transcend unity. 
With (|T71) and ( IT9l) . we get an upper bound for Aa, 

Aa < cMexp(-/3i^M), ^ = l(a - 7) > 0, (20) 

where c is some constant. The upper bound of Aa helps us determine an upper bound for 
the probability of finding a{t) deviating away from the mean a by a distance larger than e. 
Actually, following Reimann 29 1, using the Chebyshev inequality j32|, we have 

A^a 

Proh{\a{t) -a\> t) < (21) 
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For a fixed value of e, tlie upper bound decreases exponentially with the size of the system 
according to fl20l) . It then follows the statement above. 

Here some comments are worthy. Though in the derivation above we have in mind a 



sudden quench, it is easy to see that the conclusion actually applies to any type o: 
(e.g., the Hamiltonian can be changed continuously over some period, as in 13 



quench 



1^, or 



quenched multiple times as in Sec. |IV] below), as long as after some point the Hamiltonian is 
never changed again. The reason lies in that the Frobenius norm of the density matrix p{t) 
is conserved under unitary evolutions, and thus is independent of the historical or the final 
values of H{t), but is determined entirely by the initial state. As for the operator A, only the 
properties of semi-positive-definiteness and boundedness are used. Thus similar conclusions 
can apply to other operators such as a^a^afcOfc and a\a\aiai, or operators in other models. 
Finally, it should be mentioned that the conclusion relies on the fact that the quantity in 
Eq. (fTSj) is bounded by some exponentially decreasing function, which is the case only at 
finite temperatures < oo). At zero temperature, the quantity in Eq. (fTSl) is always equal 
to unity and thus the problem is still open. 



IV. A SECOND QUENCH: TYPICALITY 

It is shown in Fig. H] that after a finite transient time, the physical variables equilibrate 
to their average values exhibiting minimal fluctuations. Moreover, it has been proven that 
the amplitudes of the fluctuations will decrease exponentially with the size of the system. 
Therefore, the observation is that the system, described by the density matrix p(t), is almost 
indistinguishable from a system described by the time-averaged density matrix p, as far as 
the simple realistic physical variables are concerned. This is remarkable. Because though 
p(t) evolves unitarily and suffers no loss of information of pi, it behaves as if it were fully 
decoherenced. The question is then, to what extent can we hold onto this proposition? Is it 
possible to distinguish p{t) and p, or p(ti) and p(t2) (^i 7^ ^2), by some means? Motivated 
by this problem, we have considered the scenario of giving the quenched system a second 
quench. That is, after the first quench at t = which changes U from Ui to f/j^, at time 
t = ti, the system is quenched again by changing the value of U from Uf^ to Uf^, which is 
then held on forever. The concern is, would the long-time dynamics of the system depends 
on the specific time ti? 
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(■1 {Ui,.u„) .(-1,1) 




W) (ir,,M,,l .{2.1) 



(r) ('^'/,,'^'/J =(10,2) 



(J) (f„,I//,) -(0,!) 




FIG. 5: (Color online) The distance D between the matrices 0, and pt^ (upper panels) and the 
time- averaged values of {a\ak) (lower panels), as functions of the time of the second quench ti. 
The dashed lines in the lower panels indicate the average values of the corresponding solid lines, 
i.e., values given by Vt [see Eq. (125p ]. The initial state is the same as in previous Figures. The 
parameters (C//i, Uj^) are shown in the inserts. 
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(b) (t'c.-t'/J ={2.1 
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(J)("/„l';,) -(2,1) 



FIG. 6: (Color online) The figure-of-merit of recurrence ii as a function of time. Also shown are 
D{Q,,pt^) and (a|.afc)^^. Note the correlation between the three in (a) and (c). In (a) and (b), the 
initial state is the same as in previous Figures, i.e., {M, N,q, Dg) = (9,9,1,2700), Ui = 1, and 
/3i = 0.3. In (c) and (d), the initial state is of (M, N, q, Dq) = (6, 10, 1, 497), Ui = 1, and = 0.3. 
The values of {Uf-^,Uf^) are given in the inserts. 

Denote the Hamiltonian associated with f/jj Hf^. The density matrix of the system 
later is given by p{t) = e~*^^'2(*^*i^p(i(:i)e*^-^2(*~*i) (for t > ti). As before, we are interested 
in the fong-time averaged value of p(t), 

Pi, = lim i / dtp{h + t), (22) 

since it has been shown and proven above that the dynamics of the system is to a large extent 
captured by the time-averaged density matrix. Here the subscript indicates the dependence 
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(d) {Uj„Uf,] =(0,2) 
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FIG. 7: (Color online) Time evolution of the population on the k = Bloch state (agao) M\- Other 
/c's show similar behavior and thus are not shown. The figures correspond to those in Fig. [5] one-to- 
one. The initial state is the same as in Fig. [T] and Fig. SI In (a)-(d), the different ti's investigated 
are (0,0.35,20,40,60,70,80,90), while in (e)-(h), the different ti's are (0,0.35,5,10,15,20). Note 
that in each figure, the black and green lines correspond to ti = and 0.35, respectively. 



on the time ti. It is also useful to define the average of with respect to t 



1; 



Q = lim — 

T->oo T 



— lim 



(23) 



The second equality means that VL is actually the time-averaged density matrix associated 
with an initial state p [see Eqs. (jl]) and ^] and a Hamiltonian ifjj- O^^^ purpose of defining 
Vt is to set a reference state independent of ti. 

To gain an overall idea of the dependence on ti of the long-time dynamics, we have 
studied the distance between pt^ and Vt 28], and the time-averaged value of (a^Ofc), 



(4afc)t, = lim - / dt-tr{p{ti+t)alak) 
= tr{pt^a\ak), 



(24) 



as functions of ti. Note that the average value of {ci\cik)^^ with respect to ti is given by fi, 

I rT 

lim- / dti{alak)i^=tr{nalak). (25) 



This is another reason for defining VL. The quantities DiVt, pt^) and {a\ak)t^ shown in 
Fig. [5l Eight pairs of [Uf^^Uf^) are examined with the same initial condition as in Fig. [TJ 
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We see that for all cases with f//^ 7^ 0, both D and (0^^/0)4^ set down to their average values 
quickly. However, for the special case of Uj-^ = 0, both D and {a\ak)ti display repeated 
recurrences, without any sign of equilibration. The situation is the reverse of that in Fig. HJ 
where (a^Ofc) does not show any fluctuations in the case of f//^ = 0. 

This phenomenon is due to the recurrence of the density matrix p{t) to pi 33j. From 
Eq. (|3]), we see that in the representation of {l^/'m)}, the mn-th off-diagonal element of p{t) 
rotates at an angular frequency of — . In the generic case of f//^ 7^ 0, the energy 
gaps E^ — E^ are quite random and incommensurate, and thus recurrence of the density 
matrix is rare. More precisely, the span between two times when all the matrix elements of 
p{t) get (nearly) in phase again is extraordinarily large. On the contrary, in the special case 
of Uf^ = 0, all eigenvalues and hence all the energy gaps E^ — E^ are integral combinations 
of the few basic frequencies Uk-, and thus the probability of recurrence is much higher. To 
demonstrate that the sharp peaks in Figs. [5t and [Sji are due to recurrences of the density 
matrix p{t) to pi, we define the figure-of- merit of recurrence, 

Rit) = '^'"•"C:;" , (26) 

2—/m,n Pmn 

where the prime means the summation is over {m,n) such that E^ 7^ E^^. It is clear that 
< -R < 1 and R = 1 when and only when all the off-diagonal elements get in phase. 

In Figs. |6^ and |6]d, which share the same parameters as Figs. |5t and [5^ respectively. 



we have shown R(ti) together with D{Q,pt-^) and {ci\cik)t^- -'■^ -^^S- Ek, we see that every 
time D{Vt,pt^) and {a\ak)f-^ get close to their values at ti = 0, R{ti) shows a peak. In 
other words, there is a strong positive correlation between R{ti) and D{VL,pt^) and {ci\cik)t^- 
In comparison, in Fig. (6)3, R{ti) drops quickly from unity to less than 0.2 and remains 
low all the time, and in turn D{VL,pt^) and {ci\cik)t^ do not show any recurrence. To further 
consolidate the connection between the recurrence of pit) and that of DiVt, ptj and (a^a^)^^, 
we have considered the case of M = 6. In this case, if f//i = 0, all the basic frequencies Uk 
are commensurate, and thus there exist perfect recurrences, as shown in Fig. [6b. There we 
see clearly that -D(f2,pfJ and (a^afc)^^ return to their original values at ti = periodically, 
and this happens when and only when R returns to unity. However, once Uf-^ 7^ is set 
nonzero (see Fig. |6li) and thus the commensurability of the energy gaps is destroyed, the 
situation returns to that in Fig. [6]d. 
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The fact revealed in Fig. |5] and Fig. [6] is quite interesting. The long-time dynamics of 
the system is sensitive or insensitive to the exact time when the second quench is applied, 
depending on whether the intermediate Hamiltonian Hf^ is integrable {Uj^ = 0) or non- 
integrable (f//^ 7^ 0). In the integrable case, {ol.(ik)f^ exhibits large fluctuations and repeated 
recurrences. The system retains the memory of the initial state under the control of the 
Hamiltonian Hf^. By contrast, in the non-integrable case, (0^0^)^^ go over to their average 
values (predicted by Q) after a transitory period, showing little dependence on ti afterwards. 
Combined with Fig. HJ the picture is that p{t) evolving under the control of a non-integrable 
Hamiltonian, not only yields the expectation values of a^ak as if it were p, but even responds 
to the second quench as if it were p. 

In Fig. [TJ we have checked this picture by studying the real time evolution of (a\ak) with 
/c = under the double-quench scenario. The eight figures shown correspond to those in 
Fig. [5] respectively. For each pair of {Uf^,Uf^), we have studied the evolution of {a^aQj for 
several different values of ti. We see that in all the cases with Uf-^ 7^ 0, as long as ti is larger 
than the transient time, which can be roughly read from Fig. [5l the later evolution of {a^ao) 
is quantitatively independent of ti. On the contrary, in the case with Uf^ =0, the later 
values of (ajao) vary wildly for different values of ti. 

Here it is instructive to combine Fig. H] and Fig. [7] and compare. In the f//j 7^ cases. 



there is a sense of typicality 



35 



36|. The density matrix p(t) governed by Hf^ is surely 
non-stationary. However, for p{t) at different times, they yield almost the same expectation 
values for the observables, and moreover, they share almost the same response to the same 
quench. In the case of f//^ = 0, what Fig. [7] reveals is a good complement to that in Fig. HI 
It demonstrates that it is inappropriate to say that the system thermalizes in this case, 
even though the density matrices and expectation values of the observables agree — since 
according to one's everyday experience, a system in thermal equilibrium should not show 
any time dependence. 

V. CONCLUSIONS AND DISCUSSIONS 

We have studied the quench dynamics of the Bose-Hubbard model both analytically and 
numerically. The issues of thermalization and equilibration are investigated thoroughly. 
On the thermalization side, which concerns whether the quenched system behaves like 
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a canonical ensemble, it is found that this is the case only for small- amplitude quenches 
(at least for the finite-sized system investigated). However, the time-averaged density ma- 
trix does manifest many interesting features in different regimes. These features are self- 
consistently understood after a study of the overlaps between the eigenstatcs of Hi and 
Hf^. Here we would like to say that it is urgent and would be very helpful to develop some 
analytical tools so that some general relations between the eigen-systems of Hi and Hf^ can 
be established. These tools and relations would also be useful to determine whether the 
non-thermalization phenomenon observed is just a finite-size effect. 

On the equilibration side, which is about whether physical observables relax to station- 
ary values without appreciable fluctuations, the result is that this is indeed the case for 
quantities as {a\ak) which are of most interest. Moreover, it is proven analytically that 
for these quantities the fluctuations in time will decay exponentially with the size of the 
system. Therefore, the overall picture is that generally the system equilibrates but without 
thermalization . 

The second quench reveals something more subtle. First, the subsequent dynamics de- 
pends or not on the second quench time ti according to Uf^ = or not. The underline 
reason is the recurrence or not of the initial density matrix, which in turn has its root in 
the eigenvalue statistics of the Hamiltonian Hf^ . This effect leaves us the impression that a 
non-integrable Hamiltonian has more "dephasing power" than an integrable one. Possibly 
it can be a tool to check the integrability of a Hamiltonian. Second, in the case of C//^ 7^ 0, 
it is found that the system described by p{ti) responds to the second quench as if it were p 
for ti larger than the transient time. This means that we can take the equilibration more 
serious — p(ti) and p not only yield almost the same expectation values for the generic phys- 
ical variables but also yield almost the same dynamics after a quench. Moreover, the fact 
that the transient time is short indicates that the intermediate Hamiltonian Hf^, which is 
non-integrable, is effective in "dephasing" the initial density matrix. In another perspective, 
the dynamics of the system is sensitive to the fluctuations of U. This has the imphcation 
that in future experiments, accurate control of U would be a necessity to interpret the results 
correctly. 
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